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1.  Introduction 


Accurate  simulation  of  fracture  in  large-scale  applications  is  hampered  by  both  model  and 
numerical  algorithmic  deficiencies.  These  are  often  coupled.  The  sophistication  of  a  material 
model  can  be  limited  by  available  numerical  methods,  so  algorithm  improvements  are  often 
necessary  before  failure  model  improvements  can  be  implemented.  A  particular  algorithmic 
feature  explored  in  this  report  is  allowing  the  fracture  criterion  to  evolve  based  on  the  state  of 
adjacent  material. 

It  is  usually  easier  to  propagate  an  existing  crack  than  to  nucleate  a  new  one,  and  the  crack 
propagation  rate  is  generally  controlled  either  by  crack  tip  kinetic  processes  or  by  stress 
redistribution  at  the  crack  front.  Fracture  models  are  traditionally  considered  to  be  local,  which 
means  that  all  of  the  information  used  to  progress  the  material  failure  is  contained  in  the  close 
proximity  of  the  material  point.  While  this  is  accurate  if  the  crack  tip  is  spatially  resolved  and 
the  effects  of  the  stress  singularity  are  captured  at  neighboring  material  points,  in  typical 
large-scale  analyses,  the  numerical  discretization  is  comparatively  coarse  and  the  appropriate 
stress  and  strain  concentrations  are  not  generated  ahead  of  the  unresolved  cracks.  Further,  in  a 
local  model,  a  numerical  quadrature  point  is  unaware  that  fracture  may  exist  at  a  neighboring 
quadrature  point. 

The  goal  in  this  work  is  to  provide  infonnation  about  the  existence  of  fracture  to  neighboring 
elements  in  a  finite  element  code  and  to  adjust  the  failure  criterion  accordingly.  The  approach  is 
similar  to  the  algorithm  employed  in  Lagrangian  simulations  by  Wilkins1  and  more  recently  by 
Hohnquist  and  Johnson2.  The  development  will  be  in  the  context  of  an  Arbitrary  Lagrange 
Eulerian  (ALE)  finite  element  code,  although  the  algorithm  is  applicable  to  other  numerical 
techniques. 


2.  Approach 


The  approach  is  a  simple  nonlocal  scheme  in  which  elements  containing  fractured  material 
inform  their  neighboring  elements  of  the  earliest  time  when  the  fracture  may  arrive  at  the 
neighbor  locations.  This  enables  two  nonconventional  fracture  model  features:  the  failure 
criterion  can  be  changed  in  the  vicinity  of  a  crack,  and,  through  the  timing  of  the  criterion 
change,  the  fracture  propagation  rate  can  be  controlled  from  element  to  element.  While  the 


'Wilkins,  M.  L.  Mechanics  of  Penetration  and  Perforation.  Int.  J.  Engng.  Sci.  1978,  16,  793-807. 

2 

“Holmquist,  T.  J.;  Johnson,  G.  R.  A  Computational  Constitutive  Model  for  Glass  Subjected  to  Large  Strains,  High  Strain 
Rates  and  High  Pressures',  Report  No.  18.12544/023;  Southwest  Research  Institute:  San  Antonio,  TX,  2010. 
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technique  is  nonlocal,  it  is  implemented  in  an  operator-split  approach  and  is  minimally  disruptive 
to  the  finite  element  code.  The  method  has  many  components  implemented  within  the  traditional 
material  model  evaluations  and  a  nonlocal  component  that  involves  resetting  a  history  variable  in 
neighboring  elements. 


2.1  Computations  During  the  Constitutive  Model  Evaluation 

The  algorithm  is  illustrated  using  simple  maximum  principal  stress  failure  criteria.  There  are 
two  distinct  failure  conditions.  One  is  for  propagating  failure  in  elements  adjacent  to  material 
that  has  already  fractured,  Fprop,  and  the  other  is  for  nucleating  failure  in  regions  remote  from 
existing  cracks,  Fnuc.  Generally,  Fnuc  >  Fprop,  and  the  nucleation  threshold  can  be  much  higher 
than  the  propagation  threshold  for  ceramics  and  glasses.  The  active  failure  criterion  is  assigned 
by  a  Boolean  state  variable  within  each  element,  Acrit.  Failure  within  the  element  progresses 
when  the  maximum  principal  stress  from  an  elastic  stress  update,  crPmax,  exceeds  the  stress 
threshold: 

@Pmax  F  (1  f extent)-  (1) 

F  represents  the  threshold  from  whichever  failure  criterion  is  active  (Fprop  or  Fnuc).  An 
extent-of-failure  variable  f extent  is  introduced  to  degrade  the  threshold  strength  within  the 
element  and  control  the  stress  reduction.  It  is  essentially  a  damage  variable,  but  in  this  simple 
demonstration  implementation,  the  elastic  properties  are  not  degraded  with  damage. 

The  evolution  of  the  failure  extent  variable  within  the  constitutive  model  is  intended  to  represent 
crack  propagation  across  an  element.  If  the  crack  propagates  at  a  speed  Ccrk,  the  failure  extent  is 
updated  by 

dfenent  =  ^  * ,  (2) 

Llen 

where  L len  is  a  characteristic  element  length  and  dt  is  the  time  increment.  The  failure  extent 
variable  is  bounded,  0  <  fextent  <  1. 


The  failure  is  treated  as  being  isotropic,  so  there  is  not  a  unique  way  to  reduce  the  stress  such  that 
the  maximum  principal  stress  satisfies  equation  1 .  Here,  the  pressure  and  von  Mises  equivalent 
stress  at  the  time  of  initial  failure  are  saved  as  history  variables  patFaii  ar|d  cftFaib  respectively. 
If  the  pressure  is  tensile  (negative)  following  an  elastic  stress  update,  the  pressure  is  reset  to 
satisfy 


V  —  PatFailC^-  fextent)- 


(3) 


No  pressure  adjustment  is  necessary  if  the  inequality  is  satisfied  after  the  elastic  stress  update.  If 
the  pressure  is  compressive,  it  is  not  reset.  Failed,  dense,  solid  material  can  support  pressure,  so 
degrading  a  compressive  pressure  would  not  be  consistent  with  the  physics.  The  deviatoric 
components  of  the  stress  tensor  are  reduced  by  the  radial  return  algorithm,  if  necessary,  such  that 
the  von  Mises  stress  satisfies 


o 


vm 


—  ®atFail 


(1  fextent)- 


(4) 
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As  with  the  pressure,  the  von  Mises  stress  is  not  reset  if  it  already  satisfies  the  inequality 
following  the  elastic  stress  update.  With  the  exception  of  the  multiple  failure  criteria,  equations 
1-4  are  a  common  means  for  controlling  the  failure  rate  within  an  element. 

The  failure  criterion  is  set  by  Boolean  state  variable  Acrit.  This  variable  defaults  to  1  initially, 
representing  undamaged  material  with  failure  criterion  Fnuc .  It  is  initially  set  to  0  in  elements 
containing  significant  defects  that  would  fracture  at  the  lower  criterion,  and  set  to  0  during  the 
course  of  the  calculation  if  an  existing  crack  could  have  propagated  into  the  element.  This 
element-to-element  propagation,  described  in  the  following  section,  pennits  control  of  the  crack 
speed. 


2.2  Nonlocal  Aspects  of  the  Algorithm 

The  nonlocal  aspect  of  the  model  is  determining  when  to  alter  the  failure  criterion.  Each  element 
is  given  a  state  variable  tarrivah  representing  the  earliest  time  at  which  a  nearby  crack,  traveling 
at  the  rate  Ccrk,  could  arrive  at  the  element  centroid.  When  the  current  time  in  the  simulation  t 
meets  the  condition 

t>tarrival (5) 

Acrit  is  set  to  0,  triggering  the  use  of  the  Fprop  failure  criterion.  The  subtracted  tenn  in 
equation  5  accounts  for  the  crack  arriving  at  the  element  perimeter,  given  that  the  tarrivai 
variable  is  the  arrival  time  at  the  element  centroid.  The  tarrivai  variable  is  initially  set  to  a  value 
several  times  the  time  it  takes  a  crack  to  propagate  across  an  element.  If  failure  does  not 
progress  in  an  element  during  a  time  step,  the  tarrivai  variable  in  that  element  is  incremented  by 
the  time  step  size.  This  keeps  the  arrival  time  well  ahead  of  the  current  time  to  prevent 
premature  switching  of  the  failure  criterion.  Setting  of  Acrit  and  the  incrementation  of  tarrivai 
by  dt  are  performed  during  the  constitutive  model  evaluation. 


The  time  of  arrival  variable  tarrivat  is  subject  to  a  more  dramatic  reset  in  one  of  two  ways.  First, 
if  the  failure  extent  fextent  is  zero  (indicating  that  the  material  has  not  begun  to  fracture),  and  the 
failure  criterion,  equation  1,  is  just  met  for  the  first  time  during  the  current  time  step,  tarrivai  is 
set  to  the  current  simulation  time  plus  the  time  it  would  take  the  crack  to  arrive  at  the  element 
centroid: 


t 


arrival 


t-dt  +  -^2. 

2  Ccrk 


(6) 


The  adjustment  by  dt  puts  the  crack  inside  the  element  boundary.  During  this  time  step,  fextent 
is  incremented  by  equation  2  for  consistency  with  the  arrival  time.  Since  f extent  is  no  longer 
zero,  this  type  of  a  reset  can  only  occur  once  for  each  element.  This  reset  is  done  in  the  context 
of  the  constitutive  model  calculation. 
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The  second  way  tarrivai  is  reset  is  external  to  the  material  model  when  the  code  has  access  to  all 
of  the  current  element  history  variables.  For  each  element  in  which  fextent  >  0  (called  the 
“seed”  element  here),  the  centroid-to-centroid  distance  Dcentroid  is  calculated  to  all  adjacent 
elements  that  have  not  yet  begun  to  fail,  fextent  —  0-  The  earliest  time  for  arrival  of  the  crack  at 
the  centroid  of  these  undamaged  neighbor  elements  is  then  set  by 


C crk 

where  t*arrival  is  the  arrival  time  for  the  “seed”  element,  which  has  fextent  >  0-  The  minimum  in 
equation  7  is  necessary  since  multiple  failed  neighbors  can  reset  tarrivat  independently  and  the 
shortest  time  is  desired.  The  reset  is  only  performed  if  fextent  =  0,  so  the  minimum  arrival  time 
in  an  element  is  not  reset  in  this  nonlocal  computation  if  failure  has  already  begun. 

A  subtle  point  for  the  time  incrementation  of  tarrivai  is  that  its  value  is  increased  by  the  time  step 
in  the  constitutive  model  evaluation  if  failure  does  not  progress.  In  other  words,  the  dt 
increment  occurs  when  failure  does  not  progress  even  if  the  material  is  partially  failed, 
fextent  >  0-  The  purpose  of  tarrivai  in  an  element  with  an  existing  crack  is  to  monitor  the  time 
when  the  crack  would  finish  traversing  the  element  and  enter  a  neighbor.  For  a  stalled  crack,  the 
time  the  crack  would  exit  the  element  must  be  increased.  By  incrementing  tarrivai  for  stalled 
cracks,  the  arrival  time  in  all  of  the  surrounding  undamaged  elements  will  also  increase  during 
the  nonlocal  tarrivai  reset. 

2.3  Efforts  to  Mitigate  Advection  Issues 

Eulerian  or  ALE  calculations  move  material  across  element  boundaries.  The  algorithms  are 
designed  for  accurate  advection  of  smooth,  continuous  fields,  and  applying  these  same  methods 
to  the  discrete  variables  associated  with  failure  models  can  cause  problems.  Discrete  state 
variables  can  change  abruptly  across  element  boundaries.  Advection  algorithms  average  the  state 
of  material  entering  an  element  with  the  state  of  the  remaining  material,  and  the  smeared  result  is 
often  not  meaningful  for  the  algorithm.  Integer  variables  like  Acrit  must  be  either  0  or  1 .  The 
advection  result  may  not  be  proper  even  for  real  valued  history  variables,  such  as  f extent  •  This 
represents  the  fractional  extent  a  crack  has  propagated  across  an  element.  It  is  a  decimal  value  in 
the  element  containing  the  crack  tip,  and  it  must  be  0  in  elements  ahead  of  the  crack  and  1  in 
elements  behind  the  crack.  It  makes  no  sense  to  be  0.8  in  one  element  and  0.05  in  the  neighbor 
ahead  of  the  crack.  Hence,  some  additional  checks  and  calculations  are  developed  in  the  attempt 
to  improve  the  behavior  of  the  model  during  advection.  The  following  approach  is,  ultimately, 
not  satisfactory,  but  it  represents  some  progress  in  making  the  advected  fields  consistent. 

The  strategy  is  to  develop  the  midpoint  crack  arrival  time  state  variable  into  a  longer-range, 
smooth  field  that  will  advect  more  accurately.  Calculations  based  off  of  this  state  variable 
should  be  improved.  The  tarrival  variable  is  propagated  spatially  by  extending  the  seed  elements 
in  calculations  using  equation  7  to  include  all  elements  rather  than  only  elements  that  have 


t, 


arrival 


Min  (j- arrival’  f 


arrival 


4 


begun  to  fail.  The  tarrivai  variable  in  elements  adjacent  to  a  fracture  is  set  as  previously 
described,  and  this  serves  as  the  seed  for  setting  tarrivai  in  the  next-near  neighbor  elements  and 
so  on.  This  constructs  a  smooth,  long-range  gradient  field,  and  tarrivai  is  still  set  only  in 
elements  that  have  not  yet  begun  to  fracture.  That  keeps  tarrivai  in  actively  failing  elements  near 
the  current  time.  There  is  a  lag  built  into  the  propagation  algorithm  due  to  the  step-by-step 
propagation  of  information.  However,  the  incrementation  of  tarrivai  by  dt  advances  an 
established  gradient  field  uniformly  in  time,  so  the  lag  is  evident  only  where  the  gradient  is  being 
established  or  modified. 


The  value  of  tarrivai  is  capped  at  several  element-to-element  propagation  times  above  and  below 
the  current  time,  since  the  gradients  are  only  needed  in  the  vicinity  of  the  current  time.  Capping 
the  value  at  the  lower  end  also  circumvents  problems  where  bits  of  material  fractured  much 
earlier  in  the  calculation  can  contact  pristine  material  through  advection  and  significantly  lower 
the  tarrivai  value. 


In  Eulerian  or  ALE  codes,  there  is  often  the  opportunity  to  adjust  the  material  state  variables 
following  advection  to  bring  the  material  model  into  a  consistent  state.  During  this  adjustment, 
the  smooth  tarrival  values  are  used  to  redefine  the  more  step-like  f extent  value. 


u 


extent 


=  < 


0 

1 

\+(t-  tarrivai  )  ^ 

z  Llen 


tarrivai  >  1  +  \  ^  ~  € 
tarrivai  <  t  -  —  h  € 

otherwise 


(8) 


The  e  is  a  small  “fuzz”  value  intended  to  keep  f extent  values  above  the  numerical  noise.  In  the 
latter  case  of  equation  8,  f extent  is  subsequently  checked  to  ensure  that  it  lies  between  0  and  1. 
A  tolerance  8  is  placed  on  this  check  to  truncate  results  that  are  close  to  the  bounds.  The  effect 
of  this  tolerance  is  described  at  the  end  of  section  3.3. 


The  fracture  model  and  nonlocal  algorithm  setting  the  state  variables  were  implemented  in 
Lawrence  Livermore  National  Laboratory  code  ALE3D.3  The  nonlocal  method  is  quite  similar 
to  element-to-element  lighting  time  algorithms  for  explosive  bum,  so  much  of  the  infrastructure 
was  already  in  place.  It  is  anticipated  that  implementation  would  be  similarly  straightforward  in 
other  codes  with  similar  features. 


Nichols,  A.  L.,  Ed.  Users  Manual  for  ALE3D:  An  Arbitrary  Lagrange/Eulerian  2-D  and  3-D  Code  System',  Lawrence 
Livermore  National  Laboratory:  Livermore,  CA,  2009. 
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3.  Results 


The  ability  of  the  model  to  control  the  failure  propagation  rate  is  illustrated  through  plane-strain 
crack  growth  simulations  in  which  limiting  fracture  propagation  speeds  are  specified  as  2  and 
4  mm/ps.  The  geometry  of  the  brittle  target  plate  is  a  100  x  100  mm  square  with  a  20-mm  deep 
crack  on  the  right  side,  as  shown  in  figure  1 .  The  element  size  is  0.25  mm  in  the  x  andy 
directions,  in  both  the  plate  and  the  in  projectile.  The  plate  is  impacted  by  a  10-mm-wide  by 
20-mm-long  block  of  steel  at  0.1  mm/ps. 


Figure  1.  Initial  plane-strain  configuration 

showing  the  target  plate,  the  projectile, 
and  the  initial  crack  at  the  right  side. 

The  plate  has  a  density  of  3.9  cm  and  is  modeled  with  a  Murnaghan  equation  of  state  with  a 
reference  bulk  modulus  of  228  GPa.  The  pressure  dependence  of  the  bulk  modulus  is  2.0  and  the 
Gruneisen  parameter  is  set  to  2.  The  shear  modulus  is  152  GPa.  The  critical  principal  stress  for 
fracture  nucleation  is  set  to  30  GPa,  which  is  never  reached  in  these  simulations.  The  critical 
principal  stress  for  propagating  existing  fractures  is  set  to  200  MPa. 

The  compressive  stress  wave  initiating  at  the  impact  site  traverses  to  the  right,  where  it  reflects 
off  of  the  surface  as  a  tensile  wave.  The  tensile  wave  interacts  with  the  preexisting  crack,  and 
the  crack  grows.  Given  the  simple  nature  of  the  fracture  model,  the  cracks  tend  to  follow  the 
mesh,  but  some  cracks  are  driven  hard  enough  to  run  skew  to  the  mesh.  The  obvious  mesh 
dependence  is  not  important  for  current  purposes,  as  the  goal  is  to  demonstrate  control  of  the 
fracture  propagation  rate  and  the  influence  of  advection  on  the  failure. 
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3.1  Fracture  Rate  Control 

The  ability  of  the  algorithm  to  control  fracture  rate  is  illustrated  for  Lagrangian  simulations  in 
figure  2.  Figures  2a  and  2b  show  the  fracture  configuration  at  35  and  50  |us,  respectively,  after 
impact  in  simulations  with  the  fracture  velocity  limited  to  2  mm/ps.  Given  this  speed  and  the 
time  increment  between  the  frames,  the  most  a  crack  would  have  been  permitted  to  progress 
between  these  two  frames  is  30  mm,  less  than  one-third  of  the  distance  across  the  plate. 
Examining  the  crack  propagating  to  the  left  in  figures  2a  and  2b,  the  distance  progressed  appears 
to  follow  the  growth  restriction.  Figures  2c  and  2d  show  the  fracture  configuration  at  the  same 
times,  but  with  a  limiting  fracture  velocity  of  4  mm/ps.  It  is  evident  that  the  cracks  are  growing 
faster,  as  permitted. 


Figure  2.  Fracture  configurations  at  35  and  50  (as  for  plane-strain  simulations  with  failure  propagation  rates 
restricted  to  2  and  4  mm/ps. 

Note  that  controlling  the  failure  propagation  rate  also  affects  the  failure  pattern.  Overall,  the 
slower  propagation  rate  appears  to  create  more  regions  where  the  failure  pattern  is  spread 
spatially.  Also,  the  cracks  do  not  appear  to  follow  the  mesh  as  readily  when  the  propagation  rate 
is  more  restricted.  When  the  propagation  rate  is  slower,  the  stress  is  not  relieved  as  quickly  or 
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focused  as  tightly  along  the  mesh  lines  ahead  of  the  cracks.  This  provides  the  crack  an 
opportunity  to  grow  in  other  directions  rather  than  following  the  mesh. 

In  general,  the  cracks  tend  to  propagate  along  mesh  lines  because  the  failure  criterion  is  reduced 
along  the  mesh  direction  earlier  than  it  is  reduced  along  the  diagonal. 

3.2  Arrival  Time  Field 

The  crack  arrival  time  field  is  depicted  in  figure  3  for  the  4-mm/|u,s  propagation  rate  calculation. 
The  time  is  35  ps,  which  corresponds  to  the  color  cyan  in  the  plots.  The  zoomed  region  in  figure 
3b  is  taken  from  the  upper  left  of  figure  3  a.  The  plots  show  that  the  arrival  time  field  is  smooth 
and  continuous  for  future  times.  In  the  narrow  areas  where  the  crack  is  actively  extending,  the 
color  scale  is  continuous  to  times  earlier  than  the  present.  However,  most  of  the  dark  blue 
elements  are  adjacent  to  cyan.  This  color  scale  discontinuity  indicates  crack  surfaces  that  have 
not  propagated  for  the  last  several  time  steps.  Since  the  arrival  time  is  updated  only  with  the 
gradient  algorithm  in  elements  with  no  damage,  times  in  the  past  do  not  get  set  by  the  gradient 
algorithm.  They  are  only  truncated  to  keep  the  time  from  lagging  too  far  behind  the  current  time. 


Figure  3.  Plot  of  the  earliest  possible  time  of  arrival  of  the  fracture  at  35  ps  for  the  4-mm/|is  fracture  propagation 
speed:  (a)  full  field  and  (b)  enlargement  of  active  branch  near  the  upper  left. 


While  it  would  be  beneficial  to  have  arrival  times  continuous  into  the  past  to  provide  smoother 
fields  for  advection,  attempts  to  set  the  previous  arrival  times  in  the  gradient  algorithm  created 
feedback  loops  that  altered  the  future  times  as  well.  Further  development  in  this  area  is  needed. 

3.3  Advection  Effects 

The  same  model  geometry  is  used  to  evaluate  the  effects  of  advection  on  the  solution.  For  the 
applications  and  brittle  materials  targeted  with  this  model  construct,  it  is  anticipated  that  the 
deformation  and  material  motion  through  the  mesh  will  be  small  prior  to  fracture  and  that  large 
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deformation  is  possible  after  initial  fracture.  Fracture  patterns  from  three  simulation  approaches 
are  shown  in  figure  4:  the  first  is  a  pure  Lagrangian  calculation,  the  second  a  Eulerian 
calculation  where  the  penetrator  is  moving  into  an  initially  stationary  target,  and  the  third  an 
Eulerian  calculation  where  the  target  moves  to  the  left  into  a  stationary  rod.  The  difference 
between  the  two  Eulerian  calculations  is  the  range  of  motion  of  the  material  through  the  mesh. 
When  the  target  is  initially  stationary,  much  of  the  fracture  occurs  before  the  target  begins 
moving  appreciably  to  the  right.  When  the  target  is  initially  moving,  the  cracks  all  traverse 
several  elements  of  the  Eulerian  mesh. 


Figure  4.  Numerical  diffusion  of  fracture  field  as  a  function  of  extent  of  advection:  (a)  Lagrange  simulation,  no 
advection;  (b)  Eulerian  with  initially  stationary  target,  little  advection;  and  (c)  Eulerian  with  target 
moving  into  projectile,  significant  advection. 

Fine  cracks,  many  an  element  in  width,  span  the  Lagrangian  target  in  figure  4a.  The  cracks 
maintain  their  width  and  location  in  the  target  over  time,  which  is  consistent  with  the  physics  of 
fractured  materials.  The  fracture  pattern  is  similar  in  the  Eulerian  simulation  with  the  initially 
stationary  target,  figure  4b.  The  material  is  almost  Lagrangian  because  the  target  starts  to 
accelerate  downstream  only  as  the  momentum  is  transferred  from  the  projectile.  There  is  some 
growth  in  the  width  of  the  vertical  cracks  as  the  advection  algorithm  smears  the  cracks 
horizontally.  This  advection-related  crack  diffusion  is  greater  on  the  right  side  of  the  plot  where 
cracks  formed  earlier  and  have  been  in  motion  longer.  The  cracks  also  grow  wider  as  the 
simulation  progresses. 

The  initially  moving  Eulerian  target,  figure  4c,  shows  significant  crack  diffusion  with  mesh 
motion.  The  algorithm  setting  the  minimum  time  of  arrival  in  a  smooth  gradient  did  not  prevent 
spreading  of  the  discontinuous  crack  field.  Close  inspection  of  the  results  shows  that  the  arrival 
time  algorithm  was  partially  successful.  On  the  crack  flank  that  is  advancing  into  the  mesh  (left 
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side  of  the  vertical  cracks  in  figure  4c),  the  arrival  time  gradient  moves  with  the  material  and  the 
f extent  calculated  from  equation  8  behaves  as  desired.  The  extent,  representing  the  location  of 
the  crack  flank,  is  between  0  and  1  only  in  one  element  at  a  time,  so  that  the  leading  edge  of  the 
crack  flank  is  tracked  as  a  sharp  interface.  The  error  occurs  on  the  trailing  crack  flank.  The 
minimum  in  the  algorithm  setting  the  crack  arrival  time  continually  resets  the  timing  of  the 
gradient,  and  the  gradient  stays  approximately  fixed  to  the  stationary  mesh.  Hence,  the 
downstream  crack  flank  stays  fixed  in  space  due  to  the  calculations  in  equation  8.  This  is  no 
better  than  applying  the  advection  directly  to  integers. 

If  the  trailing  crack  flank  could  be  differentiated  from  the  leading  flank,  the  algorithm  could  be 
modified  to  treat  the  upward  gradient  in  failure  times  differently,  and  the  trailing  edge  of  the 
crack  would  also  move  through  the  mesh.  While  this  infonnation  exists  for  upwinding  of 
advected  variables,  it  is  not  directly  available  in  the  data  constructs  where  the  nonlocal 
calculations  are  perfonned. 

If  the  tolerance  8,  associated  with  the  last  of  equation  8,  is  set  large  compared  to  the  advection  of 
f extent  through  the  mesh,  the  leading  edge  of  the  crack  can  also  be  fixed  to  the  mesh  rather  than 
moving  with  the  material.  Therefore,  this  value  should  be  kept  small  enough  to  remove  only  the 
numerical  round-off  from  the  computation. 


4.  Summary  and  Conclusions 


A  nonlocal  algorithm  has  been  described  for  controlling  failure  propagation  rates  in  Lagrangian 
finite  element  calculations.  It  is  an  implementation  of  a  procedure  used  by  Wilkins1  for  ceramics 
several  decades  ago  and  more  recently  by  Hohnquist  and  Johnson."  When  implemented  in  a 
massively  parallel  finite  element  code  for  Lagrangian  simulations,  the  method  is  successful  in 
controlling  fracture  propagation  rates  with  minimal  computational  overhead. 

A  method  was  explored  to  improve  advection  of  fracture  in  Eulerian  and  ALE  codes  by 
propagating  the  earliest  time  of  arrival  of  the  crack  several  elements  away  from  the  crack 
surface.  This  created  a  smooth  gradient  field  that  would  advect  more  accurately  than  spatially 
discontinuous  state  variables  associated  with  fracture.  This  technique  was  marginally  successful, 
and  it  was  detennined  that  it  is  necessary  to  treat  the  leading  edge  and  trailing  edge  of  this 
gradient  field  differently.  The  leading  flank  of  a  propagating  crack  advanced  as  intended,  while 
the  trailing  edge  did  not.  These  results  suggest  a  path  forward,  but  further  research  into 
advection  of  discontinuous  fields  is  needed.  The  current  method  is  anticipated  to  be  well  suited 
for  materials  that  become  finely  comminuted  behind  the  failure  front,  and  accurate  representation 
of  discrete  cracks  is  not  required. 
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